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Abstract 

We present a method based on a time domain simulation of wave propaga- 
tion that allows studying the optical response of a broad range of dielectric 
photonic structures. This method is particularly suitable for dealing with 
complex biological structures. One of the main features of the proposed ap- 
proach is the simple and intuitive way of defining the setup and the photonic 
structure to be simulated, which can be done by feeding the simulation with a 
digital image of the structure. We also develop a set of techniques to process 
the behavior of the evolving waves within the simulation. These techniques 
include a direction filter, that permits decoupling of waves travelling simul- 
taneously in different directions, a dynamic differential absorber, to cancel 
the waves reflected at the edges of the simulation space, a multi-frequency 
excitation scheme based on a filter that allows decoupling waves of differ- 
ent wavelengths travelling simultaneously, and a near-to-far-field approach 
to evaluate the resulting wavefield outside the simulation domain. We vali- 
date the code and, as an example, apply it to the complex structure found 
in a microorganism called Diachea leucopoda, which exhibits a multicolor 
iridescent appearance. 
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nanobiology. 



1. Introduction 

The study of tlie optical response of structures with typical sizes of the or- 
der of the optical wavelengths has gain great interest in recent years. Emerg- 
ing technologies had resulted from the study of photonic materials, which 
consist of a regular distribution of particles within a host matrix. Depending 
on the size of the inclusions relative to the operating wavelength, photonic 
materials can be classified in metamaterials (characteristic size significantly 
smaller than the wavelength) and photonic crystals (characteristic size of the 
order of the wavelength). Then, metamaterials operating within the visible 
range are nano-structured dielectric or metallic materials, which are designed 
to control the effective electric permittivity and magnetic permeability in or- 
der to obtain a specific response. In particular, the possibility of generating a 
material with negative refraction index [ij led to a great variety of promising 
applications, such as superlensing and optical cloaking 0,31,1,8. On 
the other hand, all-dielectric photonic crystals exhibit interesting properties 
that arise from the resonant scattering generated by the specific modulation 
of the refraction index in a dielectric transparent material 0, S, loj . Photons 
in this kind of medium play the same role as electrons in the atomic crystal 



lattice [id, lUj, and then the photonic system also exhibits band gaps, i.e., 
forbidden bands for the propagation of waves. These artificial materials can 
be designed for specific purposes such as optical switches, Bragg filters or 



photonic crystal fibers [12|, |13|, [14 



Another growing research field based on dielectric photonic structures is 
the study of natural structural color, which is responsible for the iridescent 



appearance that exhibit a broad diversity of animals and plants |15l . Il6l . |17 
Structural color is produced by the selective reflection of light incident on 
the microscopic structures present in the cover tissues of biological organ- 
isms. Optical mechanisms such as interference, diffraction and scattering are 
involved to achieve colorful patterns or metallic appearance. These colors 
usually appear considerably brighter than those of pigments, although they 
often result from completely transparent materials [l8|, |l9|. Unlike artificial 
photonic materials, the geometry and distribution of these natural media is 
usually extremely complex, and the simulation of their electromagnetic re- 
sponse require versatile and accurate tools. The study of this phenomenon 
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contributes to understand different behavioural functions of living species 
such as thermoregulation and camouflage and, at the same time, inspires 
new developments of artificial devices. 

A large variety of rigorous electromagnetic methods to obtain the op- 
tical response of a given photonic structure are available in the literature. 
Among these methods, we can mention the modal method 



coordinate-transformation methods [25|, l26| and the integral method [27|, |28|, 



29] . These approaches are very efficient for the accurate determination of the 
optical response of corrugated interfaces and periodic gratings of canonical 
shapes. However, in most cases this kind of methods are not suitable for 
dealing with highly complex structures. 

Another way of studying the electromagnetic response of complex nano- 
structures is by means of computer simulations. A very spread approach is 
the Finite Difference Time Domain method (FDTD) introduced by Taflove 



et al. around 40 years ago [30|. This method is based on the Yee algorithm 



3l| and consists in numerically solving six coupled vector equations obtained 



from Maxwell's equations in the time domain. The FDTD is a very powerful 
method and has been improved during the last decades to account for a great 
variety of problems in electrodynamics. However, it results heavily time- 
consuming, and requires large computer resources and even parallelization 



for very large simulation spaces [32 



In this paper, we present a very simple simulation method that allows 
studying the propagation of electromagnetic waves in a dielectric medium of 
arbitrary refraction index distribution. Due to its simplicity, the evolution of 
the propagating waves can be easily visualized on a conventional computer 
during runtime. One of the highlights of the proposed method is its versatility 
to obtain the optical response of an arbitrary dielectric photonic structure, 
which can be introduced within the simulation by means of a diagram or a 
photograph in a form of a digital image or bitmap. Then, the structure to 
be studied can be easily generated by means of any photo-editor software. 
In the case of natural photonic structures, an electron microscope image of 
the real specimen can be used. 

In Sec. 2 we summarize the basic concepts of the simulation. Within 
the framework of the developed method, in Sections 3-6 we present a set 
of techniques which permit us to control and analyze the behavior of the 
waves within the simulation. We implemented a direction filter that per- 
mits decoupling waves travelling simultaneously in different directions, and 
also allows determining the field of energy flux in any type of wave (Sec. 
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3). In Sec. 4 we present an active system to cancel waves reflected at the 
edges of the simulation space, which allows simulating boundary conditions 
that represent an unbounded virtual space. This technique is independent 
of the wavelength and also works for any waveform, and then, it is capable 
of handling a multi-frequency excitation. Unlike other techniques to simu- 
late perfectly absorbing boundary conditions, the proposed approach does 
not need to establish a flnite thickness layer, and this constitutes one of its 
main advantages since it avoids wasting simulation space and time for that 
purpose. A multi-frequency excitation scheme based on a fllter that allows 
decoupling of waves of several wavelengths travelling simultaneously is also 
shown in Sec. 5. Such a system increases the computing speed since it 
avoids making a sequential frequency sweep to obtain the spectral response 
of a given structure. Finally, in Sec. 6 we show how to obtain the far fleld 
(outside the simulation space) from the near fleld. As an example, in Sec. 7 
we apply the simulation method, including the whole set of developed tech- 
niques, to obtain the optical response of the photonic structure present in 
the tissue of a microorganism. Concluding remarks are given in Section 8. 



2. Description of the simulation 

The simulation developed in this work is based on the method presented 



m 



33| and reproduces the propagation of transverse mechanical waves. The 



correspondence between photonics and the behavior of mechanical waves have 



been already reported in the literature [3J|. The physical model consists in 
a two-dimensional array of p x q particles of mass m contained in the x — y 
plane. Each particle is joined to its four neighbors by means of elastic springs 
of elastic constant k and separated by a distance d. The movement of the 
particles is constrained to the z-axis, which is normal to the plane of the 
two-dimensional array, as shown in Fig. [H The net force on each particle 
is null when it is located at z = 0. A wave (that will propagate along the 
X — y plane) is generated by applying an external force along the 2;-axis to 
certain particles. The dynamic evolution of the system along time is obtained 
by means of an iterative algorithm based on an integration method of flnite 



33 



timesteps 

For large p and q, the array of particles can be considered as a continuous 
medium representing an elastic membrane subjected to a tension per unit 
length T = F^/It, where Fe is the elastic force and It is the one-dimensional 
transversal section. The simulated membrane has a surface mass density 
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Figure 1: Particle array representing the physical model of the simulation. 

fi = niNm/s, where m stands for the mass of the mesh element and Nm is 
the number of particles contained in a region of area s. For a medium with 
ground density mass yUo, the speed of the waves is given by vq = \fTjJiQ. 
Therefore, in a region with a mass density fi > fio the speed of the waves is 
V < Vq- Making an analogy with optics, regions with mass density no can 
be identified with vacuum, i.e., a medium of refraction index uq = 1, while 
a region with an arbitrary mass density fi corresponds to a medium with 
a real part of the refraction index n = a/ /i//io- This approach permits to 
reproduce optical phenomena involving dielectric materials illuminated by 
transverse electric (TE) polarized light in two-dimensional configurations, 
and this has been verified for different systems and incidence configurations. 

One of the remarkable features of this method is the possibility of defining 
the simulation domain by means of digital images or bitmaps. Each pixel in 
the image represents the position of the particle in the array. In this manner, 
a digital image of px q pixels automatically defines the size of the simulation 
domain. An appropriate constant cXp given in [nm/pixel] links the size of an 
object in the digital image, in pixels, with the actual physical size of the 
sample to be simulated, which is measured in nanometers. 

Three bitmaps of equal size are defined within the simulation code to 
introduce different characteristics of the structure and the illumination con- 
ditions, which are denoted as M (mass density), D (damping) and E (ex- 
citation). The grey levels in the M bitmap represent the mass distribution 
within the array, which should be assigned using an adequate linear conver- 
sion function. Bitmap D encodes the damping constant of each point of the 
array, which allows introducing an attenuation constant in the medium, and 
bitmap E is introduced to specify the excitation, i.e., the particles on which 
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the external force is applied. Following the analogy with the propagation 
of electromagnetic waves in an optical medium, the M bitmap determines 
the refraction index distribution in space, bitmap D specifies the regions 
where there is absorption, and bitmap E specifies the location and shape of 
light sources. The grey level matrices M, D and E are related to the matri- 
ces Mphys, Dphys and Ep^ys, respectively, which contain the values of mass, 
damping and excitation measured in physical units, by: 

Mphys = mo + Mrup, (1) 

Dphys = D^p, (2) 

Ephys = rp{E- 128), (3) 

where mo is the minimum mass value, which represents the refraction index 
of vacuum. The proportionality constant rup in ([1]) has units of [kg/gl] and 
the constant /ip has units of [N sec/(m * gl)]. Notice that [gl] is the grey 
level unit within a scale of 0-255 grey levels in which the digital images are 
represented. Vp in eq. ([2]) is a constant of units of [N/gl], which transforms 
the value of grey level provided by the bitmap E to a. value of force. In 
eq. ([3]), a value E = 128 indicates that no force is applied on the particle. 
Then, values over 128 are interpreted as positive forces and grey levels values 
under 128 are interpreted as negative forces ^]. The harmonic excitation is 
introduced as 

Et = Ephys sin {u Tn + , (4) 

where Et is the applied external force, u is the angular frequency of the ex- 
citation, if is the initial phase and is a time adapting constant in units 
of [sec/cycle] that converts the integer number of iteration cycles n^. -which 
define the timestep-, into a physical time variable. The product r„ Uc repre- 
sents the discretized time variable, and the product Ud = ojTn in eq. (jl]) can 
be regarded as a digitized angular frequency measured in [rad/cycle]. 



The Nyquist-Shannon sampling theorem [35| states that the frequency of 
the signal sampling must be at least twice the highest frequency component 
of the signal, in order to preserve the alternating nature of the external 
excitation after the sampling. In our case, the sinusoidal waveform of the 
externally applied force has a period 27r, and then it should be sampled as 
minimum at twice its frequency, that is, every tt radians or less every iteration 
cycle. This implies that should be smaller than vr rad/cycle. On the other 
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hand, the adapting constants r„ and CTp are related by 

Tn = CTp, (5j 

where Vd is the digitized speed of the waves (measured in [pixels/cycle]) and 
Vphys is the physical speed of the waves measured in [nm/sec]. In the case 
of optical waves in vacuum, Vphys = vq corresponds to c = 2.99792 x 10^^ 
nm/sec. Therefore, the digitized angular frequency Ud can be expressed as: 

ujd = '^—^- (6) 

The above expression implies that once the optical frequency u (or the optical 
wavelength A) is fixed, the Nyquist-Shannon criterion requires 

VdCTp < — . (7) 
OJ 



On the other hand, the Courant-Friedrichs-Lewy condition [36| imposes 
that Vdo, the digital counterpart of the maximum allowed speed vq, must 
satisfy Vdo < 1 pixel/cycle. That is, the maximum allowed digitized wave 
speed that guarantees the stability of the simulation is one pixel per cycle of 
iteration. A speed beyond this value would cause the simulation to diverge. 
In other words, the dynamical information can be transferred to a maximum 
distance of one pixel, i.e., from one pixel to the next one, in each iteration 
cycle. This limit is not related to the computer processor speed or the physics 
of the modelled system. Instead, it ensures an internal logical consistency of 
the numerical method. 

The proposed method allows studying the response of any two-dimensional 
distribution of refraction index in a very easy and practical manner. The re- 
fraction index distribution can be artificially generated using any available 
computational design tool for digital image edition, or it can be obtained 
from a digitized electron microscope image of a real physical structure. Re- 
cently, Kolle et al. implemented an interface for the MEEP package (an 
FDTD implementation) js^], which permits introducing the profile of the 
diffracting structure via binary images based on refraction index contrasts in 
SEM or TEM images 

Unlike other numerical methods in which the mathematical expression of 
the problem is known and the computer is used as a numerical integrator 
to obtain the solution, the idea behind the proposed simulation method is 
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the use of the computer as a generator of a virtual environment where the 
physical differential law is used to make the system evolve along the time, 
as it would evolve in the real physical world. The final solution is not pre- 
determined by any mathematical expressions, but by the configuration and 
causality of the natural evolution of the introduced physical law. This way of 
thinking the problem has the advantage of naturally reproducing all the phe- 
nomena derived from the characteristics of the simulated medium and from 
the physical law involved. In this case, the simulated physical system is able 
to reproduce the whole family of linear wave phenomena, such as refraction, 
diffraction and interference, without having explicitly introduced any wave 



equation and/or boundary conditions [33 



3. Direction filter 

As mentioned above, the proposed simulation has the advantage of gen- 
uinely reproducing the phenomena that would appear in a real physical sys- 
tem. However, at the same time, this advantage is the source of the main 
drawback of the method: its uncontroUability. For instance, spurious re- 
flected waves will naturally appear at the edges of the simulation space. 
Incident and reflected waves will be superimposed, and this superposition 
makes it impossible to distinguish the amount of energy carried out by each 
one of the waves. The energy density then results from the interference of the 
incident and reflected waves. Besides, as it occurs in any experimental setup, 
within the simulation is not possible to determine the normal reflectance of a 
given sample under normally incident light, because the detector cannot be 
placed in the path of the incident light beam. And any artificial perturba- 
tion introduced in the movement of the particles within the mesh in order to 
decouple both waves, would affect the original (and actual) response of the 
system. 

Fortunately, in contrast to what happens in the physical world, all the 
dynamical information of the system is known at every time in the simulation. 
This allows extracting this information in order to be processed in different 
ways. In the general case, multiple waves will be travelling simultaneously in 
different directions, producing a complex movement, and one would like to 
decouple these different components from the complex collective movement 
of the particles within the mesh. A method used to isolate a wave travelling 
in a given direction should be based on a natural way of processing the 
information in order to be able to work automatically for any wavelength. 
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amplitude, direction, waveform and, also, in the presence of an unknown 
number of other waves travelling in different directions. 

At a first sight, this does not seem a trivial task and therefore the essential 
concept that represents a wave must be identified. A travelling wave of any 
wavelength, amplitude, direction or waveform has the following property: 
within a small region of space (small enough to be considered homogeneous), 
the waveform does not change in time, and only its phase changes at a known 
rate given by its phase velocity. In this context, we present a very simple 
and natural direction filter {DF) that allows extracting a wave travelling 
in a given direction and with a given sense. The filter is based on a wave 
subtraction technique and only uses as input parameter the speed of the wave 
to be filtered. 

In order to present the direction filter in a clear way, in the following sub- 
section we describe its formulation for one-dimensional problems and then, 
we present its extension to more dimensions. 

3.1. Mathematical formulation of the DF in ID 

Suppose we have certain time-evolving function A{x, t). Now, we propose 
the following operator: 

T^+\A{x,t)] = A{x,t + At) - A{x-vAt,t), (8) 

where w is a constant representing a speed. If and only if A{x, t) is a wave 
travelling at speed v towards the +x direction, it must satisfy 

A{x,t + At) = A{x-vAt,t). (9) 

Substituting ([9]) in ([8]) we obtain 

J^^^\A{x,t)]=A{x-vAt,t)-A{x-vAt,t)=0 Vt. (10) 

J-"*^"*"^ is called the positive DF operator acting on the wavefield A(x,t) that 
cancels waves travelling towards the +x direction. Correspondingly, J-'^~'^ = 
A{x, t + At) — A{x + vAt, t) is the negative DF operator, and it cancels waves 
travelling towards the —x direction. In other words, the DF operator cancels 
a wave A{x, t) travelling in a given direction by subtracting from it the same 
wave but evaluated in a previous instant and in a position displaced by an 
amount Ax = vAt from its present position. This is schematically shown in 
Fig. [3 
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Figure 2: Geometric representation of the action performed by the positive DF operator. 
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Let us now evaluate the effect of the positive DF operator in a more gen- 
eral case (the analysis for the negative DF operator is completely analogous). 
Suppose that A{x, t) = B~^{x, t) + B~{x, t) is the superposition of two waves 
of arbitrary shapes travelling with speed v towards opposite directions +x 
and —X, respectively. Taking into account that 

A{x - vAt, t) = - f At, t) + B-{x- vAt, t) (11) 

and 

A{x, t + At) = B^{x, t + At) + B-{x, t + At), (12) 
by applying to this new function, we get: 

J^^+\A{x,t)] = B'^{x,t+At)+B-{x,t+At)-B+{x-vAt,t)-B-{x-vAt,t). 

(13) 

Since B~^{x,t) and B {x,t) are waves travelling towards the +x and —x 
directions, respectively, they satisfy 

B^{x,t + At) = B^{x -vAt,t) (14) 

and 

B-{x,t + At) = B-{x + vAt,t). (15) 
Therefore, the application of the positive DF operator results to be 

J^W[A(x,t)] = B~{x + vAt,t) - B-{x-vAt,t). (16) 

As expected, B~^{x,t) is completely cancelled. By calling x' = x + vAt, eq. 
(|T6|l can be rewritten as 

T'-^^lAix' -vAt,t)] = B-{x,t) - B-{x,t-2At), (17) 

which represents a wave travelling towards the —x direction. By dividing 
both sides of eq. ffTTl) by 2 At and taking the limit At — )• 0, eq. ffT7|) can be 
expressed as: 

-r-l'4-^r./ M dB^ix + vAtA) 

^W[A(a;,t)] ^ 2 At ^—^^ ^. (18) 

Expression f|T8l) reveals the effect of the DF operator on a wave of general 
shape travelling in a non-filtered direction. For small At -compared with 
the time period of the higher harmonic component of B^ {x, t)-, the filtered 
wave is proportional to the time derivative of the original wave. 
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3.2. The Direction filter in more dimensions 

In general electromagnetic scattering problems we have two- or three- 
dimensional wavefields, and then in this subsection we generalize the DF to 
more dimensions. 

Consider a scalar wavefield A„(r, t) in a space of n dimensions. The 
general expression for the DF operator ([8]) in M" is 

J-^^'-^) [A„(r, t)] = A„(r, t + At) - A„(r - v5, t), (19) 

This operator filters waves travelling with phase speed |v| in the direction 
of V, with a characteristic time delay At = 5. Note that bold letters repre- 
sent vectors. To find out the effect produced by the DF operator on waves 
travelling in directions different from the filtering direction defined by v, 
we consider a wavefield A2{v,t) in M^, in which case the DF operator f ll9p 
becomes 

J-^^''^)[A2(r,t)] = A2(r,t + At)-A2(r-v5,t). (20) 
If the wavefield is a two-dimensional plane wave 

A^iv^t) = Ae'^-^''-^-*\ (21) 

where k^„ is the wave- vector, v^^, = u;/|k^| is the velocity of the wave and 
uj its angular frequency, substitution of f l2T|) into fl20|) yields 

\J^!^'^\a)\ = ^|e-^'^^-e'''"'^=°'^(")|, (22) 

where 1^2^'^'' (a) | is the complex amplitude of the filtered wave and a is the 
angle between the propagation direction (k^^,) and the filtering direction (v) 
(see Fig. [3D. 

To quantify the performance of the DF, we define the relative attenuation 
/ia of the wave as 

^M) = 1 - |J•i'''')(a)|/|J•^■')r^^ (23) 

where | J^^^^ '^'' |™^^ stands for the maximum value of I J'a^^'^^ (a) | for a G [0, 360°] . 

Figure H] shows the relative attenuation of a plane wave as a function of 
a. As expected, the maximum attenuation is obtained for a = 0, that is, 
when the propagation direction of the plane wave coincides with the filtering 
direction. Conversely, the attenuation is minimum for the wave travelling 
in the opposite direction, that is, for a = 180°. In Fig. |5] we show the 
performance of the DF for a circular wavefront, which can be regarded as a 
superposition of plane waves propagating along all directions. The amplitude 
of the obtained filtered wave is shown in grey levels. 
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Totally filtered wave 





a [Degrees] 

Figure 4: Attenuation of the plane wave as a function of the angle a between its direction 
of propagation and the filtering direction. 




Figure 5: Angular effect of the DF on a circular wavefront. 
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3.3. Direction filter implementation within the computer simulation 

As presented above, the DF is characterized by two parameters: v, which 
determines the speed and direction of the wave to be filtered, and 5, the 
characteristic delay. As in any digital simulation, the space-time domain is 
discretized and then, this discretization imposes certain constrains on the 
DF parameters. 

From eq. fl2U]) it becomes evident that the filtering operation will be 
effective if the spatially displaced A2{r — ^r5,t) has exactly the same shape 
as A2(r, t + At). Since in the discretized domain the minimum separation 
distance is 1 pixel, and the simulation allows a maximum digitized wave speed 
VdQ < 1 pixel/cycle (according to the Courant-Friedrichs-Lewy condition), 
the travelling wave advances a distance = vt = Vdo x ArZc < 1 pixel in one 
iteration cycle (Aric = 1), and therefore, A2{r,t + At) will have the same 
shape as A2{y — w5,t) exactly after 1/vdo iteration cycles, with (l/v^o) £ ^■ 
During intermediate iteration steps, the wave will take interpolated values 
between A2{y — \5,t) and 742(r,t + At) that will not exactly match the shape 
of neither of them. 

According to this, the necessary requirement for a good performance of 
the DF within the simulation is that the speed of the waves must be set 
to Vdo = 1 pixel/Anc, with Anc G Z being the integer number of iteration 
cycles in which the wave advances a distance of exactly 1 pixel. Therefore, 
the characteristic delay of the filter must be 5 = Ari^ = 1 pixel/f^o for 
a spatial shifting of A2{y — vS,t) equal to 1 pixel. In our simulations we 
usually set the second allowed digitized speed Vdo = 0.5 pixels/cycle, (i.e. 
Aric = 2). Therefore, in this case 6 = 2 cycles for a spatial shifting of 1 
pixel in A2{r — v5, t). Although we could eventually use spatial shiftings 
larger than 1 pixel, larger shiftings lead to larger values of the characteristic 
delay 6, during which the shapes of A2{r — v5, t) and A2{r,t + At) would 
be more affected by the numerical errors, thus decreasing the quality of the 
directional filtering. 

It is worthwhile to mention that the above conditions are valid for the im- 
plementation of the DF in the orthogonal directions x and y, which coincide 
with the rows and columns of the discretized array of particles comprising 
the simulation medium, as shown in Fig. [1] In the present paper we use four 
DFs corresponding to the directions +x, —x, +y and —y. From a practical 
programming point of view, the DF is applied to the waves evolving in the 
two-dimensional array defined by the simulation, called main plane. How- 
ever, the resulting directionally filtered waves are stored in secondary planes. 
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Figure 6: (a) Source emitting gaussian waves in directions +x and —x. (b) Secondary 
plane of the DF that cancels waves traveUing towards the +x direction, (c) Secondary 
plane of the direction filter that cancels waves travelling towards the —x direction. Blue 
arrows show the propagation direction of the wavefronts. 

These secondary planes show time-evolving waves that are filtered images 
of the waves evolving in the main plane, and therefore, the DF does not 
affect the physics of the system. In this case, we have four secondary planes 
that store the filtered waves travelling in the +x, —x, +y and —y directions, 
respectively. 

To illustrate the behavior of the DF, in Fig. EJ^a) we show the main 
plane of the simulation of a linear source emitting gaussian waves in both 
directions ±x. Figures |6t^b) and[6]^c) show the secondary planes that contain 
the wavefield of Fig. [6](a) after the application of the DF that cancels waves 
travelling towards the +x and the —x direction, respectively. 

3.4- Determination of the energy flux vector field 

One of the main applications of the DF presented in the previous subsec- 
tions is the method for evaluating the energy flux of any simulated wavefield. 
This method makes use of the waves filtered in the orthogonal directions ±x 
and ±y. Let us call Ir^^, I^-, Iy+ and /y_ the intensity (calculated as the time 
integration of the square of the amplitude) of the waves travelling towards 
the +x, —X, +y and —y directions, respectively. Then, the quantities 

(p^ = h+ - 4_ (24) 

and 

0y = h+ - h~ (25) 
are proportional to the energy flux along the x and y directions, respectively. 
Consequently, the vector field F defined as 

F = (0.,0,) (26) 
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is proportional to the energy flux vector field. 

According to (ITS]) , the filtered intensities are not equal to the intensity 
of the original wavefield, but proportional to its derivative. However, to 
determine the direction of the energy flux -and not its absolute value-, it is 
enough to compute the relative intensities between the filtered waves. 

As an example, let us calculate (px for a wave Z[r,t) travelling along the 
+y or the —y direction. Since 

J-t''\z{r,t)]=J-t''\z{r,t)], (27) 

i.e., the DFs for the +x and —x directions produce the same attenuation 
on waves that propagate along a direction orthogonal to their corresponding 
directions of filtering (see eq. ( 122|) ). according to (!24|) (j)^ becomes 

0. = 4+ - = J i^t'^^ydt - J i4-'''^ydt = 0. (28) 

As expected, there is no energy flux along the x direction for the evaluated 
wave. 

It should be mentioned that the DF developed here could also be ap- 
plied to experimental data. For instance, the method of Fourier transform 



profilometry presented in [39[ permits digitalizing the evolution of surface 
water waves along time. Therefore, the technique proposed in this paper 
also enables the calculation of the energy flux field for real systems. 

4. Dynamic differential absorber 

It is well known that any computer wave simulation can only reproduce 
the propagation of waves in a finite domain. Therefore, the waves arriving to 



the edges of the simulation domain will be naturally reflected back [33|. In 
order to simulate open boundaries, it is necessary to artificially cancel all the 
waves reflected at the edges of the simulation space. For this purpose, several 
methods have been proposed 33, 4^, 41 1, each of them having its advantages 
and disadvantages. One possibility to avoid these reflections is to place a 
slab of an absorbing medium adjacent to the edges of the domain. In its 
basic implementation, this approach reduces the reflected waves, but it does 



not completely cancel them |33|. Besides, the absorbing region produces an 
unnecessary increment of the size of the simulation space and, consequently, 
the computation time also increases. 
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A highly improved version of the procedure described above is the so 



called perfectly matched layer (PML) method [30|, |4l|]. This method is the 
most widely spread technique to cancel reflected waves, and it consists in in- 
troducing an absorbing anisotropic layer at the edges of the simulation space. 
Within this layer, the differential wave equation is modified by including a 
special transformation that produces a rapid attenuation of the wave as it 
propagates. Although this method results to be very efficient, it requires a 
layer of finite thickness to allow the decay of the waves. Besides, waves trav- 
elling parallel to the layer are not attenuated by the PML method, producing 
the accumulation of non-realistic energy in that region. 

In this Section we present an alternative approach to simulate the open 
space by cancelling reflected waves at the edges of the simulation domain. 
The method, called dynamic differential absorber (DDA), is based on an in- 
tuitive concept, and its main advantage is that it does not require a layer 
of a given thickness to cancel the waves, i.e., it produces the absorption of 
the wave within a layer of infinitesimal thickness, and this saves computation 
space and time. On top of that, the method automatically cancels waves of 
any amplitude, shape or frequency with the same efficiency, and this consti- 
tutes a great advantage that enables the use of multi-frequency excitation, as 
shown in Section 5. As in the case of the DF, only the velocity of the incom- 
ing waves should be specified. It is important to mention that this method 
can be applied to waves of any nature, i.e., mechanical, electromagnetic, and 
potentially even to quantum waves governed by Schrodinger's equation. 

In the following subsection we present the basic idea by means of the 
ID formulation of the DDA, called here the simple dynamic differential ab- 
sorber (SDDA). Then, we develop the adaptive dynamic differential absorber 
(ADDA), which is a generalization to two or more dimensions, and in which 
case the direction of the incoming waves must be taken into account. 

4.I. Simple dynamic differential absorber (SDDA) 

Consider a wave propagating along the x-axis towards the +x direction, 
and let us examine the movement of this wave at two fixed points x = Xa and 
shown with the blue dots in Fig. [71 As the wave passes through, 
the amplitude At, at Xb will take the value of Aa at Xa, but with a time delay 
given by At = Ad/v^, where Ad = Xb — Xa and is the phase velocity of 
the wave. Within this context, the term "amplitude" refers to the magnitude 
of the wavefield at a given point. 
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Figure 7: Amplitude of a propagating wave at two fixed points located in x = Xa and 

X — Xb- 

The key idea is the following: if we artificially control the amplitude Aj, of 
the wave at xi, and make it move in such a way that Af, copies the amplitude 
Aa with the right time delay, the behavior of the propagating wave for x < 
will not be affected as it passes through. In fact, if point located at 

the edge of the simulation domain, the incoming wave will not be affected at 
this point and it would continue propagating as if there were no boundary. 
Therefore, the SDDA consists in controlling the movement of a point located 
at Xb, called the absorbing point, according to the behavior of the incoming 
wave at Xa, called the reading point. Mathematically, we can exprees this 
procedure as 

4*^) =4*0), (29) 

where A^^^ is the amplitude of the wave at Xb in a certain instant ti, Aa°^ 
is the amplitude of the wave at Xa in a previous instant to = — where 
At = 6d Act. In this context, 6d = l/v^ is the differential delay which 
depends on the medium characteristics via the phase velocity of the waves 
within the propagating medium. 

Although the above equations indicate that the reading point can be 
located arbitrarily close to the absorbing point, within the simulation method 
the discretized nature of the space-time domain must be taken into account, 
as already explained in subsection 13. 3l for the DF. Since the minimum distance 
between the reading and the absorbing points is of one pixel, we also set the 
allowed digitized wave speed Vdo = 0.5 pixels/cycle (see subsection I3.3p . and 
then 6d = 1/fdo = 2 iteration cycles. If, for instance, a ID simulation space 
is M pixels long, the rightmost pixel m = M is forced to move as: 

A^M^ = (30) 
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where Aj^^ is the amphtude of the wave at the pixel M — 1 stored two 
iteration cycles before the present cycle ric- Similarly, the leftmost pixel 
m' = 1 is forced to move as: 

where A^2^''~'^^ is the amplitude of the second pixel stored two iteration cycles 
before the present cycle Uc- 

4-2. Adaptive dynamic differential absorber (ADDA) 

If the simulation domain is two-dimensional, the SDDA presented in the 
previous subsection would only be effective for waves propagating normally 
to the edge of the simulation space (or whose wavefronts are parallel to 
the edges of the domain). In order to develop a direction-sensitive method 
that could properly absorb waves propagating with different directions, small 
corrections must be introduced into the differential delay used to cancel the 
incoming waves. 

The effective wavelength Aefr of the wave that arrives at the edge of the 
simulation space is given by 

Aeflf = A^/ cos(a), (32) 

where A^ is the actual wavelength of the incoming wave and a is the angle 
between the direction of propagation and the normal to the edge of the 
simulation space. Then, the effective phase velocity Ves -the phase velocity 
of the wave along the direction normal to the edge of the simulation space-, 
depends on the angle of incidence and is given by: 

VeE = J— = ^csf, (33) 

with u = 27r/, k^s being the effective wavenumber and / the frequency of 
the wave. This is schematically shown in Fig. |H1 Consequently, the effective 
differential delay 6ef[ becomes 

1 cos(a) 

des = — = = ddCos{a), (34) 

which implies that 6cs is always smaller than 6d- 

To determine 6es, the reading point must provide information not only 
about the time- varying amplitude of the wave to be cancelled, but also about 
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Figure 8: Effective wavelength Acff and effective phase velocity Wcff of an obliquely incident 
wave forming an angle a with the normal to the simulation space edge. 

its direction of propagation, given by a. This angle can be obtained by 
evaluating the energy flux vector at the reading point, using the method 
described in subsection 13.41 For example, for the right hand-side edge of the 
simulation space, a is given by 



where (f)x and (f)y are the x— and y— components of the energy flux, defined 
in eqs. and fl2^ . In practice, to clearly establish the direction of the 
incoming wave before it is affected by the edge of the simulation space, the 
information about the direction of the incoming waves is taken from points 
located between 2 and 4 pixels away from the edge of the simulation space, 
where the absorbing points are located. On the other hand, the discretized 
nature of the simulation must be taken into account. If the wave speed is set 
to Vdo = 0.5 pixels/cycle, 5eff is a real number between 2 and for a between 
and 90°, respectively. Since the number of cycles 6es must be an integer, 
the resulting value At must be approximated to the closest integer value. 
This imphes that the allowed discretized delays {6cs = 2, 1 and 0) will only 
match the required delay for three angles of incidence, that in this case are 
a = 0°, 60° and 90°. Therefore, in order to increase the number of allowed 




(35) 
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Figure 9: Delay vs. angle of incidence: discretized (solid line) and calculated (dashed 
line) . 

discretized delays, the wave speed can be decreased, and set, for instance, 
to Vdo = 0.1 pixels/cycle. In this case, 5eff will be bounded between 10 and 
0, for a between 0° and 90°, respectively. Then, the number of discretized 
delays is increased, and the error between the calculated and the allowed 
integer delays is minimized, as shown in Fig. |9l Increasing the number of 
discretized delays improves the effectiveness of the ADDA. Then, to increase 
the number of discrete delays we can either increase the distance between the 
reading and the absorbing point, or decrease the wave speed Vdo- Depending 
on the requirements of speed and usable space size of the simulation, the best 
choice should be made. Figure [TO] shows the location of the absorbing points 
(in blue) and the reading points (in red) for a two-dimensional space. The 
distance Ad (magnified for clarity) between the reading and the absorbing 
points is also indicated. 

Despite the difference between the calculated and the discrete values of 6d 
introduced for certain incidence angles by the discretization process, it was 
verified that these errors do not significantly affect the performance of the 
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Absorbing points 
Reading points 

Figure 10: Location of the absorbing and the reading points for a 2D simulation space. 
The distance Ad is also indicated. 

ADDA, even for 6d = 2 (second maximum digitized wave speed allowed) and 
Ad = 1 (maximum available space), as will be shown in the next subsection. 
Making an analogy with electronics, the proposed method behaves like an 
active device, which reads an input signal and reacts accordingly to return 
a post-processed output signal. In this sense, this method is different from 
most available absorbing methods, which behave like passive devices whose 
response does not take into account the input signal. 

4-3. Validation 

In order to evaluate the performance of the ADDA, the reflectance, de- 
fined as the ratio of the reflected to the incident intensity, was calculated for 
a gaussian beam incident on the edges of the simulation space with different 
angles. 

In all cases, the distance between the reading and the absorbing points was 
set to Ad = 1 pixel. Figure [TT] shows the reflectance as a function of the angle 
of incidence for Sd = 2 (blue line), (5d = 4 (green line) and Sd = S (red line). 



23 



0.15 




0.05 



40 

a [DEG] 



80 



Figure 11: Numerical experiment with the ADDA: Reflectance vs. angle of incidence a 
for 5(i = 2 (blue line), = 4 (green line) and (5c; = 8 (red line). 



It can be observed that as 8^ is increased, the overall reflectance decreases, 
that is, a better performance of the absorber is obtained. For = 2, there 
are peaks of relatively high reflectance for 39° and 69°. In the case of (5^ = 4 
the reflectance also shows peaks, but less intense than for = 2. These 
peaks are located at the angles of incidence for which the difference between 
the calculated and the available integer delays is maximized. As the number 
of integer delays is increased (by increasing 6d), a better matching between 
the calculated and the available integer delays is obtained (see Fig. IH]), and 
the intensity of the reflectance peaks is gradually reduced, as observed for 
5d = 8m Fig. EH 

Figure [12] shows the attenuation produced by the ADDA expressed in 
decibels (dB), and calculated as x = —10 \ogiQ{R), where R is the reflectance. 
It can be noticed that for the angles for which the matching between the 
integer and the calculated delays is better, the attenuation values reach up 
to 44.2 dB for 6d = 2, 64.5 dB for 5^ = 4 and 49.4 dB for 6d = 8. As a 
reference, the maximum attenuation values obtained with the PML applied 
within the FDTD method lie between 20 and 39 dB for small incidence angles, 



depending on the setting parameters of this method |41 
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Figure 12: Attenuation vs. angle of incidence a for (5^ = 2 (blue line), 5^ = 4 (green line) 
and 8d — '^ (red line). 

It can also be observed in Fig. [12] that for angles of incidence greater 
than approximately 25°, the best (larger) attenuation is obtained for 5^ = 8, 
whereas for small angles of incidence the attenuation obtained for 5^ = 4 
is better. This result can be explained by taking into account that smaller 
angles of incidence require higher values of the effective differential delay (5eff , 
as determined by In this case, there are more iteration cycles between 
the reading and the absorbing instants for 6^ = 8 than for 6^ = 4. Then, if the 
sampling resolution of the shortest wavelength involved is low, a wave that 
propagates a fraction of pixel will be represented in the next timestep by a 
natural interpolation produced by the simulation. When the wave advances 
one pixel, the accumulated numerical errors produced by interpolation will 
be higher in the case of 5^ = 8 than in the case of 6d = 4, producing a higher 
mismatch between the read waveform and the waveform that actually reaches 
the absorbing point. This produces a better performance for 5^ = 4 for small 
angles of incidence. Therefore, in order to reduce the interpolation errors 
and improve the absorbing characteristics for small incidence angles, -i.e., 
for high differential delays 5^,- the sampling resolution of the wavelengths 
comprising the wave to be cancelled must be increased, and this is done by 
decreasing the constant ap described in Section |2l 
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Figure 13: Intensity diagram of a gaussian beam forming an angle a — 20° with the 
normal to the lower horizontal edge of the simulation space for the case 6d ~ 4. (a) 
without ADDA; (b) with ADDA. 

The effect of the ADDA is graphically shown in Fig. [13] for an incident 
gaussian beam forming an angle of 20° with the normal to the edge of the 
simulation space. Figure [TST a) shows the intensity diagram without the 
ADDA, where the reflected wave can be observed. On the other hand, Fig. 
[T^ b) shows that the reflected beam is almost completely cancelled when the 
absorber is activated. 

5. Tuning filter 

The optical response of a structure is usually described by its reflectance 
and transmittance as a function of the wavelength. Therefore, in order to ob- 
tain its optical response using a simulation, the program should be executed 
for many frequency values and this could take considerable computing time, 
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especially if the desired spectral resolution is high. To avoid this problem, 
in this Section we show that a multi-frequency excitation (MFE) scheme can 
be implemented within the simulation in a quite straightforward manner. 

Within this framework, the single frequency excitation (SFE) described 
in dlj) is replaced by 



where /tot is the total number of frequencies, Ei is the excitation bitmap for 
the i-th frequency cuj and is the i-th initial phase. In what follows we set 
Ei = E, meaning a uniform frequency spectrum. 

Each point in the simulation space has a complex oscillating movement 
resulting from a mixture of the whole set of waves of different frequencies 
having unknown amplitudes, and we are interested in extracting each fre- 
quency component from the multi-frequency wavefield. One possible way to 
do this is to store the time evolution of the whole set of pixels in the simula- 
tion space in order to perform the Fourier transform for each pixel. Then, the 
spectral amplitude for each pixel can be obtained from the Fourier spectrum 
evaluated at any single frequency. The maximum number of iteration cycles 
used in each particular case depends on the size of the bitmap that defines 
the simulation space and on the setting of the digitized speed of the wave 
Vdo- In order to guarantee that the waves could reach every point within the 
simulation space, we established the following criterion: the total number 
of iteration cycles should be such that it ensures that a wave can travel 
twice the longest straight distance within the simulation space. For a bitmap 
of p X q pixels and a wave speed Vdo, the minimum number of cycles is: 



Therefore, to evaluate the Fast Fourier Transform (FFT) for each pixel, it 
would be necessary to store the complete evolution of the waves in a 3D 
matrix of size {p,q,n^c^). Hence, for a 3D simulation space, a 4D matrix 
should be stored. For example, a 2D simulation space of 1 Mpixel would 
require to store a 21 GBytes matrix (with double precision), which exceeds 
the RAM memory of today conventional computers. On the other hand, 
in the FFT algorithm the frequency domain size equals the digitized time 



/tot 




(36) 



i=l 
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domain size. Then, for a typical resolution of 10 nm in the wavelength, the 
set of points corresponding to the whole set of frequencies results to be very 
localized in the Fourier space, as shown in Fig. [TH Under these conditions, 
and due to the digitized nature of the Fourier space, the estimation of each 
spectral amplitude is highly inaccurate. In order to reduce the error, the 
size of the frequency domain should be enlarged by increasing even more the 
size of the digitized time domain n*°* given by eq. fl57|) . Figure IT^ a) shows 
the amplitude (in arbitrary units) of the time evolution of a pixel under 
multi-frequency excitation for 41 simultaneous waves of optical frequencies 
contained within the range 380 - 780 nm and using = 10000 cycles, 
Tn = 3.335 * 10^^'^ sec/cycle and the same excitation amplitudes for all the 
frequencies. The amplitudes of the detected frequencies are very irregular, 
showing that the FFT is not useful to adequately isolate the amplitude of 
each frequency under these excitation conditions and digital resolution, as 
shown in the example of Fig. [T^ b). 

To solve this problem, we introduce the tuning filter {TF), that allows ex- 
tracting a single- frequency time- varying wave from the multi-frequency wave- 
field. This filter acts as a temporal mask which is applied on the simulation 
space. After the application of the TF to the time-varying multi-frequency 
wavefield, a dynamic single-frequency wave can be visualized, allowing its 
processing to determine physical magnitudes such as reflectance or transmit- 
tance for a given frequency. 

5.1. Implementation 

The TF is based on the tuning properties of the forced damped oscillator, 
governed by the well known inhomogeneous differential equation 

--n^ + l-n + ^do^. (38) 



m dt"^ dt 

where Fext is the applied external force, z and t are the position and time 
variables, respectively, 7 is the damping constant, and Udo = ^/kdo/fndo is the 
natural frequency of the damped oscillator {kdo is the spring elastic constant 
and rrido is its mass). The TF uses the characteristic frequency response of the 
forced damped oscillator, which acts as a narrow band filter and maximizes 
the signal at resonance, i.e., when the frequency of Fext equals the natural 
frequency Udo of the oscillator. Under the condition 7 <^ Udo, the width of 
the resonance peak is given by 

Au ^ 7, (39) 
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Figure 14: Multi-frequency excitation corresponding to 41 optical wavelengths contained 
within the range 380 - 780 nm, for n^°* = 10000 cycles, t„ = 3.335 * 10"!^ sec/cycle and 
equal magnitude of the external force applied to all the waves of different frequencies: (a) 
time evolution of a single pixel; (b) the absolute value of its FFT. 
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Figure 15: Schematic diagram of the tuning filter implementation. 



the resonant frequency is given by 



(40) 



and the sharpness of the resonance peak is determined by the quahty factor 
Q, defined as 

(41) 
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Larger Q values correspond to sharper resonance peaks. 

The implementation of the TF is similar to that of the DF. The simulation 
space where the multi-frequency wavefield is evolving is considered as the 
main plane. Then, a damped oscillator is associated to each point of the main 
plane, forming a two-dimensional array of oscillators, called the secondary 
plane. As schematized in Fig. [151 the wavefield amplitude at each point in 
the main plane is used to generate a proportional force that is applied to the 
particle of mass m^o in the corresponding oscillator. The displacement of each 
particle in the secondary plane is associated to the amplitude of the filtered 
wavefield. The oscillators are tuned at the frequency to be isolated, and then, 
the value of Udo is selected accordingly. Therefore, as many secondary planes 
as frequencies to be extracted will be needed. 

The main advantage of the MFE is that the computing time is signifi- 
cantly reduced compared with the SEE. In principle, the MEE employs the 
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same time to complete the simulation of the wavefield as that used by the SFE 
for just a single frequency. For a spectrum containing /tot discrete frequen- 
cies, the total number of iteration cycles is reduced /tot times with respect 
to the number of iteration cycles required to process the same signal when a 
sequential swept of /tot single frequencies is carried out. However, the advan- 
tage in speed is compensated by a requirement of more memory to store the 
whole set of secondary planes and their respective auxiliary variables, which 
also introduces a slight delay in the computing time of each iteration cycle 
of the MFE compared with that of the SFE. This delay increases with the 
number of simultaneous frequencies explored in the MFE. More precisely, the 
computing times for a single iteration cycle within each excitation scheme 
are related by: 

^MFE = ^SFE (1 + ^MFE /tot), (42) 

where Smfe is the additional fraction of tspE required by each frequency. 
In the case of the algorithm implemented in this work, Sufe ~ 0.03879. 
In general, (5mfe could vary according to the dynamic memory allocation 
efficiency of the implemented algorithm. 

Taking into account fl42l) . to run a whole set of frequencies the MFE 
requires less computing time than the SFE if 

n%'>{l + 6MFEftot)n%, (43) 

where n*"^ is the total number of iteration cycles for the MFE and n^°g = 
/tot ■ n^°^ is the total number of iteration cycles for the SFE. 

While n^"g is only determined by the size of the simulation space and 
by /tot, in the case of the MFE n*'^ is also conditioned by the required se- 
lectivity of the TF, which depends on Au and on /tot- According to ([31]), 
the separation between adjacent discrete frequencies should be Aco > 7 to 
perform an adequate isolation of frequencies. Besides, the TF should work 
under the stationary regime, i.e., when the non-tuned frequency oscillations 
have vanished. Since the time employed by the oscillator to reach the sta- 
tionary regime is inversely proportional to 7, the time required by the TF 
to isolate the selected frequency basically depends on the desired resolution 
Au. Consequently, the value of 7 should be adjusted according to the speed 
and precision requirements of the simulation in each particular case, and this 
determines the value of n^"^- According to the above considerations, one 
could decide which scheme is more appropriate. 
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5.2. Discretization requirements and validation 

In the case of the optical spectrum (A G [380 nm, 780 nm]) discretized 
in steps of 10 nm (41 frequency values), the value of n^°^ required for an 
appropriate separation of frequencies results to be large since, in this case, 
Au = 7 is small and then the required quality factor is very high (see eq. 
f HTj) ). Moreover, to minimize the noise produced by the non-tuned compo- 
nents in the tuned frequency wave, it is even convenient to choose 7 ^ Au, 
and this implies a high value of n^^M- 

The total number of cycles n^°^ required by the TF to extract a frequency 
could, in principle, be reduced by increasing the time adapting constant r„ 
(see eq. (jlj)), i.e., by increasing Ud- As mentioned in Section [21 the Nyquist- 
Shannon criterion imposes Ud < tt. However, it was found that for high values 
oiud the TF resonates at a frequency higher than the tuned frequency. Then, 
to ensure that the TF selects the desired frequency, Ud must satisfy Ud <^ tt, 
and therefore: 

Vd (Tp = < — , 44 

CO UJ 

for ujd and u being the maximum digitized and physical frequency contained 
in the analyzed spectrum, respectively. 

In order to quantify the performance of the TF independently from other 
characteristics of the simulation, we define the relative error (in percentage) 
of the obtained intensity as a function of the explored wavelength as: 

..(A.) = 100 , (45) 

where /mfe = (^m^'^)^ ^^'^ -^sfe = (^4™^"")^ are the intensities obtained with 
the MFE for the wavelength Ae and the SFE, respectively, and and 
^max maximum amplitudes detected at the final stage of the time 

evaluation (in which the TF is in the stationary regime) for the MFE and 
the SFE cases, respectively. 

The exploration of the optical spectrum with a resolution of 10 nm is an 
extreme situation for the MFE scheme because, as mentioned above, for 41 
simultaneous frequencies a very high quality factor is needed, which leads to 
a very large number of iteration cycles n^°M- This limitation can be overcame 
by exploring the same amount of frequencies in several MFE stages having a 
larger separation between frequencies, in such a way that different frequencies 
are covered in each stage. In this manner, the required quality factor of the 
TF is greatly decreased and therefore, n**^/ is also considerably reduced. For 
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example, the optical spectrum can be explored in five MPE stages with a 
resolution of 50 nm (/tot = 9). In this case, an error of Er = 5.25% can be 
achieved with nji'^j = 50000 cycles per stage, 7 = 10""^ rad/cycle, ap = 19 
nm/pixel and Vdo = 0.5 pixels/cycle. Therefore, the total number of cycles 
required to explore the 41 frequencies will be 5n^"j^ = 250000. For 5mfe ~ 
0.03879, the MFE scheme becomes more suitable than the SEE in terms of 
computing time to explore the 41 optical frequencies, if the simulation space 
becomes larger than 1454 x 1454 pixels. This example shows that for large 
simulation spaces, the MEE (and the use of the TF) represents an advantage 
over the SEE, that can be implemented in any conventional computer. The 
decision on which scheme (MEE or SEE) is more suitable will in general 
depend on the maximum allowed error in each particular case. 

If we are just interested in the visualization of the field associated to each 
frequency, less precision is required. In this case, the whole set of frequencies 
can be analyzed in a single stage, with a higher value of 7, which considerably 
reduces the required number of iteration cycles. 

As an example, the field scattered by an opaque cylinder of diameter 620 
nm was simulated for an incident optical multi-frequency plane wave (A e 
[380 nm,780 nm], AA = 10 nm) with 7 = 10~^ rad/cycle, ap = 20 nm/pixel, 
Vdo = 0.5 pixels/cycle, and a simulation space of 150 x 150 pixels. Eigure 
[roT a) shows the resulting multi-frequency intensity diagram for the scattered 
wavefield and Eigs. [T^ b)-(d). are the intensity diagrams for several extracted 
components of wavelengths 780 nm, 570 nm and 380 nm, respectively, for 
^cM — 600. The white dashed line denotes the cylinder position. As can 
be observed, even for a very low number of iteration cycles the TF allows 
a clear visualization of the wavefield of each frequency. This result can be 
explained by taking into account that if 7 is relatively high, the non-tuned 
frequency components are rapidly damped and the filter would immediately 
be oscillating at its tuning frequency, although the stationary regime has not 
been strictly reached. 

The possibility of getting a clear and rapid visualization of the selected 
frequency components by means of the TF could also be useful to decouple 
the different frequencies that coexist in a multi-frequency field of mechanical 
waves obtained experimentally. A few thousands of frames of an experimental 
evolving multi-frequency wavefield could be enough to decouple and visualize 
the evolution of a single frequency component. This feature also makes the 
TF a valuable tool for visualization and signal processing of experimental 
data. 
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6. Near-to-far-field transformation method 



The method proposed in this work provides the near field distribution of 
a wave interacting with a given object. However, in many practical cases 
one is interested in the far-field response of an illuminated structure. On 
the other hand, there are many analytical methods that only provide the far 
field response of canonical problems, that could serve for validation and for 
appropriately setting the parameters of our simulation method. Therefore, 
in this Section we explain how the far field is obtained from the near field. 

The transformation method used to obtain the far field from the near 



field is based on Green's Theorem |30l. |42| . The remarkable advantage of this 



method is that it avoids the enlargement of the simulation domain to obtain 
the far-field information. 

If the wavefield is represented by a complex phasor Z{r,t) = Zii{v,t) + 



i Zj{r,t), according to Green's theorem 30| for a fixed time we have: 

Z,(r) = / [G(r|r') n', ■ VZ,{r') - Z,{v') n', ■ V'G(r|r')] dC\ (46) 

JCa 

where Ztiv) is Z(r, t) evaluated at a fixed t, r' is the position of a source point 
over an arbitrary contour Ca enclosing the scatterer, ria is the outward unit 
normal to the contour C^, r is an observation point outside C^, and G(r|r') is 
the Green's function, which in two dimensions is given by the Hankel function 

G{vW)='-H^i\k\v-v'\), (47) 

with i being the imaginary unit and k = 2'k/\ the wavenumber. 

To apply eq. fH6|) and calculate the far field, we need to have the complex 
near field ^t(r) for each pixel. However, the present simulation method 
provides a real scalar wavefield, and then, its imaginary part should be found. 
Taking into account that 

Zr(t) = Ar(t)e-^'^* = A-r{t) cos{ut) - i Ar{t) sin(wt), (48) 

where Ar{t) is the amplitude of the phasor at a fixed position r and u is 
the angular frequency in a steady state (when Ar{t) is constant), it is easily 
verified that 

Z,(r,t) = |{Z«(r,t)}. (49) 

Since the simulation provides Zji{r,t), we use fHOj) to calculate Zj{r,t) 
and therefore, to build the phasor for each pixel that is introduced into eq. 
fH^ to calculate the far field. 
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7. Application example 



In this Section we show an application of the proposed simulation and 
of the whole set of techniques developed in this work to obtain the optical 
response of a highly complex structure of biological origin. In particular, 
we evaluate the optical response of the peridium -a transparent protective 
layer that encloses the mass of spores- of the Diachea leucopoda (Physarales 
order, Myxomycetes class), which is a microorganism that has a characteristic 
pointillistic iridescent appearance (see Fig. [TTT a)). It has been recently 



demonstrated that this appearance is of structural origin [18|, |43[ . We use a 
transmission electron microscope (TEM) image of the peridium cross section 
of the Diachea leucopoda, shown in Fig. [TTT b). to set the M bitmap, which 
defines de refraction index distribution by means of a linear conversion of the 
grey levels of the negative image of Fig. [TTT b). The grey level (black) is 
associated to the lowest value of the refraction index (equal to unity). The 
average refraction index of the peridium was set to 1.78, which corresponds to 
a common value found in biological tissues. The optical source is a gaussian 
beam of width 2 /xm. 41 optical frequencies were explored in the range 380- 
780 nm with a spacing of 10 nm. The simulation space was set to be of 280 x 
230 pixels and cTp = 15 pixel/nm. The width of the sample obtained from the 
TEM image is 1.68 /xm and its mean thickness is 550 nm. In order to avoid 
errors produced by the finite size of the peridium image, we extended the 
biological slab at both sides by planar homogeneous slabs whose thickness is 
the average thickness of the actual image and whose refraction index is its 
average refraction index. 

Figure [TS] shows the resulting near field distribution of reflected intensity 
produced by the peridium cross section for a wavelength of 380 nm, and Fig. 
[T9] shows the total far field reflectance as a function of the observation angle 
a and of the wavelength A (the color scale bar is the same for Figs. [18] and 

m. 



8. Conclusion 

In this paper we presented an electromagnetic wave simulation method 
to obtain the time- varying fields for multi-frequency excitation and arbitrary 
shapes and refraction index distribution of dielectric objects. This method 
results particularly suitable for investigating the electromagnetic response of 
photonic structures of biological origin, that usually present a high degree 
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Figure 18: Reflected near-field intensity diagram produced by tlie peridium of the Diachea 
leucopoda obtained with the simulation method, for a wavelength of 380 nm. 
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Figure 19: Far-field reflectance (R) as a function of the observation angle a and of the 
wavelength A for the peridium of the Diachea leucopoda. 
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of complexity in their geometry as well as in the materials involved. To im- 
prove the performance of the method from the point of view of computing 
time optimization and space saving, we implemented a set of techniques that 
permit us controlling and analyzing the propagating waves within the sim- 
ulation. These methods include a direction filter, that permits decoupling 
waves travelling in different directions, a dynamic absorber, that prevents the 
reflection of waves at the edges of the simulation space in order to reproduce 
unbounded spaces, and a tuning filter that allows multi-frequency excitation. 
We also adapted a near-to-far field method to calculate the far field with a 
minimum use of computation time and allocation space. Also, possible ap- 
plications of these techniques have been suggested to analyze time-varying 
experimental wavefields. As an application example, we calculated the re- 
flectance of the transparent cover layer of a microorganism that exhibits 
iridescence, for multiple wavelengths and observation angles. 

In its present form, the proposed simulation method computes the elec- 
tromagnetic response in the case of transverse electric (TE, electric field 
perpendicular to the plane of incidence) polarized incident light. In order 
to fully simulate the electromagnetic response of complex structures with 
traslational invariance, also the transverse magnetic polarization mode (TM, 
magnetic field perpendicular to the plane of incidence) should be included. 
This would require an extension of the method to a fully vectorial formu- 
lation, that is already under development. Another aspect in which we are 
planning to work in the near future is the extension of the simulation to deal 
with three-dimensional objects. As stated above, most biological structures 
are highly complex and require a 3D model to properly account for their 
electromagnetic properties. The development of such a tool would constitute 
a valuable contribution for the study of natural photonic structures. 
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